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Caracterisation et approximation des ensembles 
de valeurs propres des matrices symetriques 

intervalles 

Resume : We consider the eigenvalue problem for the case where the input 

matrix is symmetric and its entries perturb in some given intervals. We present 
a characterization of some of the exact boundary points, which allows us to 
introduce an inner approximation algorithm, that in many case estimates exact 
bounds. To our knowledge, this is the first algorithm that is able to guaran- 
tee exactness. We illustrate our approach by several examples and numerical 
experiments. 

Mots-cles : Matrice d'intervalles, matrice symetrique, analyse par intervalles, 
valeur propre, borne sur les valeurs propres 
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1 Introduction 

Computing eigenvalues of a matrix is a basic linear algebraic task used through- 
out in mathematics, physics and computer science. Real life makes this prob- 
lem more complicated by imposing uncertainties and measurement errors on the 
matrix entries. We suppose we are given some compact intervals in which the 
matrix entries can perturb. The set of all possible real eigenvalues forms a com- 
pact set, and the question that we deal with in this paper is how to characterize 
and compute it. 

The interval eigenvalue problem has its own history. The first results are 
probably due to Deif [TU] and Rohn & Deif [M]: bounds for real and imaginary 
parts for complex eigenvalues were studied by Deif [TU], while Rohn & Deif 
[M] considered real eigenvalues. Their theorems are applicable only under an 
assumption on sign pattern invariancy of eigenvectors, which is not easy to 
verify (cf. [S]). A boundary point characterization of eigenvalue set was given 
by Rohn [30], and it was used by Hladik et al. [Ij to develop a branch & prune 
algorithm producing an arbitrarily tight approximation of the eigenvalue set. 
Another approximate method was given by Qiu et al. [55]. The related topic of 
finding verified intervals of eigenvalues for real matrices was studied in e.g. [?]. 

In this paper we focus on the case of the symmetric eigenvalue problem. 
Symmetric matrices naturally appear in many practical problems, but symmet- 
ric interval matrices are hard to deal with. This is so, mainly due to the so-called 
dependencies, that is, correlations between the matrix components. If we "for- 
get" these dependencies and solve the problem by reducing it to the previous 
case, then the results will be greatly overestimated, in general (but not the ex- 
tremal points, see Theorem [2]). From now on we consider only the symmetric 
case. 

Due to the dependencies just mentioned, the theoretical background for the 
eigenvalue problem of symmetric interval matrices is not wide enough and there 
are few practical methods. The known results are that by Deif [TD] and Hertz 
[15) . Deif [TU] gives an exact description of the eigenvalue set together with re- 
strictive its assumptions. Hertz [15] (cf. [32]) proposed a formula for computing 
two extremal points of the eigenvalue set — the largest and the smallest one. As 
the problem itself is very hard, it is not surprising conjectures on the problem 
[25] turned out to be wrong [55] . 

In the recent years, several approximation algorithms have been developed. 
An evolution strategy method by Yuan et al. |35) yields an inner approximation 
of the eigenvalue set. By means of matrix perturbation theory, Qiu et al. j27] 
proposed an algorithm for approximate bounds, and Leng & He [23] for an 
outer estimation. An outer estimation was also considered by Kolev [55], but 
for general case with nonlinear dependencies. Some initial bounds that are easy 
and fast to compute were discussed by Hladik et al. [TS] . An iterative algorithm 
for outer estimation was given by Beaumont |3]. 

This problem has a lot of applications in the field of mechanics and engi- 
neering. Let us mention for instance automobile suspension systems |28j . mass 
structures [57] > vibrating systems [TT], principal component analysis [T5], and 
robotics [J. Another applications arise from the engineering area concerning 
singular values and condition numbers. Using the well-known Jordan- Wielandt 
transformation [131 1201 125| we can simply reduce a singular value calculation to 
a symmetric eigenvalue one. 
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The rest of the paper is structured as follows: lu Sec. [5] we introduce the 
notation that we use throughout the paper. In Sec. [3] we present our main 
theoretical result, and in Sec. 2] we develop our algorithms for the problem. 
Finally, in Sec. [S] we demonstrate our approach by a number of examples and 
numerical experiments, and we conclude in Sec. [51 

2 Basic definitions and theoretical background 

Let us introduce some notions first. An interval matrix is defined as 

A := [A'A] = {Ac M™^"; A < A < ^}, 
where A, A e R™^", A<A^ are given matrices. By 

A,:^\{A + A), A^:=\{A-A) 

we denote the midpoint and the radius of A, respectively. 

By an interval linear system of equations Ax = b we mean a family of 
systems Ax = h, such that A G A, 6 G 6. In a similar way we introduce interval 
linear systems of inequalities and mixed systems of equations and inequalities. 
A vector x is a solution of Ax = b if it is a solution of Ax = h for some A £ A 
and b G b. We assume that the reader is familiar with the basics of interval 
arithmetic; for further details we refer to e.g. [3l fT4l |26] . 

Let be a family of n x n matrices. We denote the eigenvalue set of the 
family J- by 

A(J") := {A eR; Ax^ Xx, x^O, Ae T}. 
A symmetric interval matrix as defined as 

■.= {AeA \A^A^}. 

It is usually a proper subset of A. Considering the eigenvalue set A (A), it, 
generally, represents an overestimation of A(A^). That is why we focus directly 
on eigenvalue set of the symmetric counterpart, even though we must take into 
account the dependencies between the elements, in the definition of A^ . 

A real symmetric matrix A E M"^" has always n real eigenvalues, let us sort 
them in a non-increasing order 

Ai(^) > X2{A) > ■ ■ ■> A„(^). 

We extend this notation for symmetric interval matrices 

A,(A^) :={A,(A) I AgA^}. 

These sets represent n compact intervals Xi{A^) — [X^{A^), Xi{A^)], i = 
1, . . . ,n; cf. pTS]. The intervals can be disjoint, can overlap, or some of them, 
can be identical. However, what it can not happen is one interval to be a proper 
subset of another interval. The union of these intervals produces A{A^). 
Throughout the paper we use the following notation: 



RR n° 7544 



Characterizing and approximating eigenvalue sets of symmetric interval matricesb 



the ith eigenvalue of a symmetric matrix A (in a non- 
increasing order) 

the ith singular value of a matrix A (in a non-increasing 
order) 

the eigenvector associated to the ith eigenvalue of a sym- 
metric matrix A 

the spectral radius of a matrix A 
the boundary of a set <S 
the convex hull of a set <S 
the diagonal matrix with entries yi, . . . , j/rt 
the sign vector of a vector x, i.e., sgn(x) = 
(sgn(a;i), . . . ,sgn(a;„))^ 

the Euclidean vector norm, i.e., ||a;||2 = ^x'^x 

the Chebyshev (maximum) vector norm, i.e., ||a;||oo = 
max{|x|i; i = l,...,n} 

vector and matrix relations are understood component-wise 

3 Main theorem 

The following theorem is the main theoretical result of the the present paper. 
We remind the reader that the principal m x m submatrix of a given n x n 
matrix is any submatrix obtained by eliminating any n — m rows and the same 
n — m columns. 

Theorem 1. Let A G dh.{A^). There is a principal submatrix G M*^^*^ of 
such that: 

• If X = Xj{A^) for some j G {1,. . . ,n}, then 

A e {A,(Ac + diag(z)lAdiag(z)); z e {±l}^ i = l,...,k} . (1) 

• If \ = \j{A^) for some j e {1, . . . , n}, then 

Ae {Ai(Ie-diag(^)lAdiag(0)); ^ e {±l}^ i = l,...,k] . (2) 

Proof Let A € dK{A'^). Then either A = \j{A^) or A = \,j{A^), for some 
j € {1, . . . , n}. We assume the former case. The latter could be proved similarly. 

The eigenvalue A is achieved for a matrix A ^ A. It is without loss of 
generality to assume that the corresponding eigenvector x, \\x\\2 = 1, is of the 
form x = {0^,y^)^, where y € M'^ and yi ^ 0, for all 1 < i < fc, and for some 
k e {1, . . . ,n}. The symmetric interval matrix A^ can be written as 

where B'^ c M(»-fc)x("-fe), c C M("-'=)xfe, c 
From the basic equality Ax = Xx it follows that 

Cy = for some C G C, (3) 



Xi{A) 
Vi{A) 

ds 

conv<S 

diag(y) 
sgn(a;) 




x<y, A<B 
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and 

Dy — Xy for some D e D^. (4) 

We focus on the latter relation; it says that A is an eigenvalue of D. We will 
show that is the required principal submatrix and D could be written 
as in ID). 

From ^ we have that A — y'^ Dy, and hence the partial derivatives are 
dX 



ddij 



This relation strongly influences the structure of D. If j/ij/j > 0, then dij = dij. 
This is so, because otherwise by increasing dij we also increase the value of A, 
which contradicts our assumption that A lies on the upper boundary of A{A^). 
Likewise, yiyj < implies dij = d^j. This allows us to write D in the following 
more compact form 

D = Dc + diag(2)£iAdiag(z), (5) 

where z = sgn(y) G {±1}*^. Therefore, A belongs to a set as the one presented 
in the right-hand side of ([1]), which completes the proof. □ 

Note that not every X^{A^) or Xj{A^) is a boundary point of K{A^). The- 
orem [T] is also true for such Xj[A^) or Xj{A^) that are non-boundary, but 
make no multiple eigenvalue (since the corresponding eigenvector is uniquely 
determined). However, truthfulness of Theorem [1] for all Xj{A^) and Xj{A^), 
j = l,...,n, is still an open question. Moreover, full characterization of all 
Xj{A^) and Xj{A^), j = 1, . . . ,n, is lacking, too. 

As we have already mentioned, in general, the eigenvalue set of an interval 
matrix is larger than the eigenvalue set of its symmetric counterpart. This is 
true even if both the midpoint and radius matrices are symmetric (see Exam- 
ple [1]). The following theorem says that overestimation caused by the additional 
matrices is somehow limited by the intermediate area. 

Theorem 2. Let Ac.,A/\ E M"^" be symmetric matrices. Then 

conv A(A'^) = convA(A). 

Proof. The inclusion conv A(A'^) C conv A(A) follows from the definition of the 
convex hull. 

Let A £ Ahe arbitrary, A one of its real eigenvalues, and x the corresponding 
eigenvector, where ||a;||2 ~ 1. Let B :— ^{A + A^) G A^ , then the following 
holds: 

A = x^Ax < max y^ Ay = max y^ By = Xi{B) £ convA(A'^). 
I|y||2=i l|y||2=i 

Similarly, 

X = x^Ax> min y^ Ay ^ min By — Xn{B) £ comr K{A^). 
I|y||2=i l|y||2=i 

Therefore A S convA(A"^), and so convA(A) C convA(A'^), which completes 
the proof. □ 
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4 Inner approximation algorithms 

Theorem [T] naturally yields an algorithm to compute a very sharp inner approx- 
imation of K{A^)^ which could also be exact in some cases. We will present the 
algorithm in the sequel fSection l4.3p . First, we define some notions and propose 
two simple but useful methods for inner approximations. 

Any subset of S is called an inner approximation. Similarly, any set that con- 
tains S is called an outer approximation. In our case, an inner approximation of 
the eigenvalue set Xi{A^), is denoted by ^JL^{A^) — [^j,.{A^),jLi{A^)\ C \i{A^), 
and an outer approximation is denoted by a;i(A'^) — [ui^[A^ ) i{A^ )] 3 Aj(A"^), 
where 1 < i < n. 

From a practical point of view, an outer approximation is usually more use- 
ful. However, an inner approximation is also important in some applications. 
For example, it could be used to measure quality (sharpness) of an outer approx- 
imation, or it could be used to prove (Hurwitz or Schur) unstability of certain 
interval matrices, cf. [31j. 

4.1 Local improvement 

The first algorithm that we present is based on local improvement search tech- 
nique. A similar method, but for another class of symmetric interval matrices 
was proposed by Rohn [31]. The basic idea of the algorithm is to start with an 
eigenvalue, Xi{Ac), and the corresponding eigenvector, Vi{Ac), of the midpoint 
matrix, Ac, and then move to an extremal matrix in A^ according to the sign 
pattern of the eigenvector. The procedure is repeated until no improvement is 
possible. 

Algorithm [1] outputs the upper boundaries JL^{A^) of the inner approxima- 
tion [^i^A^),Jl^{A^)], where 1 < i < n. The lower boundaries, can 
be obtained similarly. The validity of the procedure follows from the fact that 
every considered matrix, A, belongs to A^ . 



Algorithm 1 (Local improvement for ii^{A^), i = 1, . . . ,n) 

1: for z = 1, . . . , n do 

2: M»(^^) = -oo; 

3: A := Ac\ 

4: while \i{A) > Ji^iA^) do 

5: £> := diag(sgn(?;j(A))); 

6: A:^ Ac + DA/^^D] 

7; M,(A^^) :=A,(^); 

8: end while 

9: end for 

10: return 'J1^{A^), i = l,...,n. 



The algorithm terminates after, at most, 2" iterations. However, usually in 
practice the number of iterations is much smaller, which makes the algorithm 
attractive for applications. Our numerical experiments (Section[5]) indicate that 
the number of iterations is rarely greater than two, even for matrices of dimen- 
sion 20. Moreover, the resulting inner approximation is quite sharp, depending 
on the width of intervals in A^ . This is not surprising as whenever the input 
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intervals are narrow enough, the algorithm produces, sometimes even after the 
first iteration, exact bounds; see fTU'. We refer the reader to Section [S] for a 
more detailed presentation of the experiments. 

4.2 Vertex enumeration 

The second method that we present is based on vertex enumeration. It consists 
of inspecting all matrices 



and continuously improving an inner approximation /J, whenever {Az) > 

7Ij(A'^), where 1 < i < n. The lower bounds, ^^{A^), could be obtained in a 
similar way using the matrices Ac — diag(2:) diag(z), where z S {±1}", and 
zi = 1. The condition zi — 1 follows from the fact that diag(z) diag(z) — 
diag(— z) Aa diag(— z), which gives us the freedom to fix one component of 
z. The number of steps that the algorithm performs is 2"^^. Therefore, this 
method is suitable only for matrices of moderate dimensions. 

The main advantages of the vertex enumeration approach are the following. 
First, it provides us with sharper inner approximation of the eigenvalue sets 
than the local improvement. Second, two of the computed bounds are exact; 
by Hertz [H] (cf. [32]) and Hertz fW\ we have that JIiiA^) = Ai(A-^) and 
li^{A^) = A„(A'^). The other bounds calculated by the vertex enumeration, 
even though it was conjectured that there were exact |29| . it turned out that 
they were not exact, in general [35] The assertion by Hertz [T6', Theorem 1] 
that ^j,^{A^) = X^{A^) and 'JI^{A^) — Xn{A^) is wrong, too; see Example |31 
Nevertheless, Theorem [T] and its proof indicate a sufficient condition: If no 
eigenvector corresponding to an eigenvalue of A^ has a zero component, then 
the vertex enumeration yields exactly the eigenvalue sets \i{A^), i = 1, . . . ,n. 

The efficient implementation of this approach is quite challenging. In or- 
der to overcome in practice the exponential complexity of the algorithm, we 
implemented a branch & bound algorithm, which is in accordance with the sug- 
gestions of Rohn [21]. However, the adopted bounds are not that tight, and the 
actual running times are usually worse than the direct vertex enumeration. That 
is why we do not consider further this variant. The direct vertex enumeration 
scheme for computing the upper bounds, 'p,^{A^), is presented in Algorithm[5J 

4.3 Submatrix vertex enumeration 

In this section we present an algorithm that is based on Theorem[TJ and it usually 
produces very tight inner approximations, even exact ones in some cases. The 
basic idea the algorithm is to enumerate all the vertices of all the principal 
submatrices of A^ . The number of steps performed with this approach is 



To overcome the obstacle of the exponential number of iterations, at least in 
practice, we notice that not all eigenvalues of the principal submatrices of the 
matrices in A^ belong to some of the eigenvalue sets Xi{A^), where 1 < i < n. 
For this we will introduce a condition for checking such an inclusion. 

RR n° 7544 
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Algorithm 2 (Vertex enumeration for i = 1, . . . ,ri) 

1: for i — 1, . . . ,n do 

2: -L^{AS) ^ \,{A,)- 

3: end for 

4: for all z e {±1}", zi = 1, do 

5: A :— Ac + Aia.g{z) A\dia,g{z)] 

6: for i — 1, . . . ,n do 

7; if A,(A) > 7Zi(A^) then 

8: := A.(A); 

9; end if 
10: end for 
11: end for 

12: return 'p.i{A^)^ i — l,...,n. 



Assume that we are given an inner approximation /Xj(A'^) and an outer ap- 
proximation u}i{A^) of the eigenvalue sets Xi{A^); that is ij,^{A^) C Xi{A^) C 
uJi{A^), where 1 < i < n. As we will see in the sequel, the quality of the output 
of our methods depends naturally on the sharpness of the outer approximation 
used. 

Let C M'^^'^ be a principal submatrix of A^ and, without loss of gener- 
ality, assume that it is situated in the right-bottom corner, i.e.. 




where C R("^'=)^("~'=) and C C r("~'=)^'=. Let A be an eigenvalue of some 
vertex matrix D £ which is of the form and let y be the corresponding 
eigenvector. If the eigenvector is not unique then A is a multiple eigenvalue and 
therefore it is a simple eigenvalue of some principal submatrix of D^] in this 
case we restrict our consideration to this submatrix. 

We want to determine whether A is equal to \p{A^) e K{A^) for some fixed 
p G {1, . . . ,n}, or if this is not possible, to improve the upper bound Ji,p{A^)\ 
the lower bound can be handled accordingly. In view of ([3]), it must be 

so that A to be an eigenvalue of some matrix in A^ . Now, we are sure that 
A e K{A^) and it remains to determine whether A also belongs to Xp{A^). 

If A < 11p{A^), then it is useless to further considering A, as it will not extend 
the inner approximation of the pth eigenvalue set. If p = 1 or A < uip_i{A^), 
then A must belong to \p{A^)^ and we can improve the inner bound 'jlp{A^) 
A. In this case the algorithm terminates early, and that is the reason we need 
(jJi{A^) to be as tight as possible, where 1 < i < n. 

If p > 1 and A > LO_p_i{A^), we proceed as follows. We pick an arbitrary 
C E C, such that Cy = 0; we refer to, e.g. [33 for details on the selection 
process. Next, we select an arbitrary B G and let 

A.^{^r g)^ (0) 
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We compute the eigenvalues of ^, and if 7ip(A"^) < Xp{A), then we set'pp{A^) :— 
Xp{A), otherwise we do nothing. 

However, it can happen that A = Xi{A^), and we do not identify it, and 
hence we do not enlarge the inner estimation iip{A^). Nevertheless, if we apply 
the method for all p = 1, . . . , n and all principal submatrices of , then we 
touch all the boundary points of A{A^). If A e dA{A^), then A is covered by 
the resulting inner approximation. In the case when A is an upper boundary 
point, we consider the maximal i £ {1, . . . ,n} such that A — Xi{A^) and then 
the zth eigenvalue of the matrix ([6]) must be equal to A. Similarly test are valid 
for a lower boundary point. 

Now we have all the ingredients at hand for the direct version of the sub- 
matrix vertex enumeration approach that is presented in Algorithm [21 which 
improves the upper bound 'JIp{A^) of an inner approximation, where the index 
p is still fixed. Let us also mention that in step |4] of Algorithm [3l the decom- 
position of A^ according to the index set J means that is a restriction of 
A^ to the rows and the columns indexed by J, is a restriction of A^ to the 
rows and the columns indexed by {1, . . . , n} \ J, and C is a restriction of A^ 
to the rows indexed by {1, . . . , n} \ J and the columns indexed by J. 



Algorithm 3 (Direct submatrix vertex enumeration for /ip(A'^)) 

compute outer approximation uJi{A^), i — 1, . . . ,n; 
call Algorithm [T] to get inner approximation iJ,^{A^), i = 1, . . . , n; 
for allJ C {1, . . . , n}, J 7^ 0, do 

decompose A^ — (^^t according to J; 

for all z e {±1}I"'I, zi = 1, do 
D -.^ Dc + diag(z) Da diag(z); 
for i = 1, . . . , I J| do 

A:=A,(i^); 
y := Vi{D); 

if A > JIpiA^) and A < ujp{A^) and e Cy then 
if p = 1 or A < uj_p_j^{A^) then 

-PpiA^) A; 
else 

find C e C such that Cy = 0; 

if Xp{A) > 7Ip(A^) then 
:= XpiA); 

end if 
end if 
end if 
end for 
end for 
end for 

return 7Ip(^'^)- 
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4.3.1 Branch &: bound improvement 

In order to tackle the exponential worst case complexity of Algorithm [31 we 
propose the following modification. Instead of inspecting all non-empty subsets 
of {1, . . . , n} in step [31 we exploit a branch & bound method, which may skip 
some useless subsets. Let a non-empty J C {!,..., n} be given. The new, 
possible improved, eigenvalue A must lie in the interval A := [pp[A^),ujp{A^)]. 
If this is the case, then the interval matrix — XI must be irregular, i.e., it 
contains a singular matrix. Moreover, the interval system 

{A'^ -XI)x^O, ||x||oo = l , 

has a solution x, where Xi — for all i ^ J. We decompose A^ — \I according to 
J, and, without loss of generality, we may assume that J — {n — \J\ + 1, . . . , n}, 
then 

AS-XI^^^'-^^ 



XI, 

The interval system becomes 

C2/ = 0, (r>^-A/)y = 0, ||2/||oo = l, (7) 

where we considered x = (O"^,?/^)^. This is a very useful necessary condition. 
If ([7]) has no solution, then we cannot improve the current inner approximation. 
We can also prune the whole branch with J as a root; that is, we will inspect no 
index sets J' C J. The strength of this condition follows from the fact that the 
system ([7]) is overconstrained, it has more equations than variables. Therefore, 
with high probability, that it has no solution, even for larger J. 

Let us make two comments about the interval system ([7]). First, this system 
has a lot of dependencies. They are caused from the multiple occurrences of A, 
and by the symmetry of . If no solver for interval systems that can handle 
dependencies is available, then we can solve ([7]) as an ordinary interval system, 
"forgetting" the dependencies. The necessary condition will be weaker, but still 
valid. This is what we did in our implementation. 

The second comment addresses the expression ||y||oo — 1- We have chosen 
the maximum norm to pertain linearity of the interval system. The expression 
could be rewritten as — 1 < y < 1 (for checking solvability of ([7]) we can use 
either normalization ||j/||oo = 1 or ||y||oo < !)• Another possibility is to write 

~1 < y < l. Hi — ^ for some i E {1, . . . , | J|}. 

This indicates that we can split the problem into solving \J\ interval systems 

Cy = 0, (D^ - XI)y = 0, -1 < y < 1, = 1 , 

where i runs, sequentially, through all the values {1, . . . , | J|}; cf. the ILS method 
proposed in [TTj. The advantage of this approach is that the overconstrained 
interval systems have (one) more equation than the original overconstrained 
system, and hence the resulting necessary condition could be stronger. Our 
numerical results discussed in Section [5l concern this variant. As a solver for 
interval systems we utilize the convex approximation approach by Beaumont 
[5]; it is sufficiently fast and produces narrow enough approximations of the 
solution set. 
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4.3.2 How to conclude for exact bounds? 

Let us summarize properties of the submatrix vertex enumeration method. On 
the one hand the worst case complexity of the algorithm is rather prohibitive, 
0(3"), but on the other, we obtain better inner approximations, and sometimes 
we get exact bounds of the eigenvalue sets. Theorem [T] and the discussion in the 
previous section allow us to recognize exact bounds. For any i e {2, . . . , n}, we 
have that if Xi{A^) < X^_-^{A^), then J1^{A^) = Xi{A^); a similar inequality 
holds for the lower bound. This is a rather theoretical recipe because we may 
not know a priori whether the assumption is satisfied. However, we can propose 
a sufficient condition: If aJj(A'^) < uj_^_j^{A^), then the assumption is obviously 
true, and we conclude JIi{A^) — Xi{A^); otherwise we cannot conclude. 

This sufficient condition is another reason why we need a sharp outer ap- 
proximation. The sharper it is, the more times we are able to conclude that the 
exact bound is achieved. 

Exploiting the condition we can also decrease running time of submatrix 
vertex enumeration. We call Algorithm [3] only for p £ {l,...,n} such that 
p = 1 or oJp{A^) < Up_-f^{A^). The resulting inner approximation may be a 
bit less tight, but the number of exact boundary points of A{A^') that we can 
identify remains the same. 

Notice that there is enough open space for developing better conditions. For 
instance, we do not know whether Jl^{A^) < ^^_^{A^) (computed by submatrix 
vertex enumeration) can serve also as a sufHcient condition for the purpose of 
determining exact bounds. 

5 Numerical experiments 

In this section we present some examples and numerical results illustrating prop- 
erties of the proposed algorithms. The experiments we performed on a PC 
Intel (R) Core 2, CPU 3 GHz, 2 GB RAM, and the source code was written in 
C++. We use GLPK v . 4 . 23 [24J for solving hnear programs, CLAPACK v . 3 . 1 . 1 
for its linear algebraic routines, and PROFIL/BIAS v . 2 . . 4 [5T] for interval 
arithmetic and basic operations. Notice, however, that routines of GLPK and 
CLAPACK[1J do not produce verified solutions; for real-life problems this may not 
be acceptable. 

Example 1. Consider the following symmetric interval matrix 

/ 1 2 [1,5]\^ 
A^ = \ 2 1 1 • 
V[l,5] 1 1 J 

Local improvement (Algorithm [T|) yields an inner approximation 

fj-i{A^) = [3.7321, 6.7843], 
fj-2{A^) = [0.0888, 0.3230], 
H^{A^) = [-4.1072, -1.0000]. 

The same result is obtained by the vertex enumeration (Algorithmic)). Therefore, 
'Jli{A^) = Xi{A^) and fi^{A^) — X^{A^). An outer approximation that is 
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needed by the submatrix vertex enumeration (Algorithmic]) is computed using 
the methods of Hladik et aL [HI [TO]. It is 

wi(A^) = [3.5230, 6.7843], 
U2iA^) = [0.0000, 1.0519], 
UsiA^) = [-4.1214, -0.2019]. 

Now, the submatrix vertex enumeration algorithm yields the inner approxima- 
tion 

fi'iiA^) = [3.7321, 6.7843], 
H'2{A^) = [0.0000, 0.3230], 
tj.'^{A^) = [-4.1072, -1.0000]. 

Since the outer approximation intervals do not overlap, we can conclude that 
this approximation is exact, that is, \i{A^) = fi[{A^), i — 1,2,3. 

This example shows two important aspects of the interval eigenvalue prob- 
lem. First, it demonstrates that the vertex enumeration does not produce exact 
bounds in general. Second, the symmetric eigenvalue set can be a proper subset 
of the unsymmetric one, i.e., A{A^) C A{A). This could be easily seen in the 
matrix 

It has three real eigenvalues 4.6458, —0.6458 and —1.0000, but the second one 
does not belong to A{A^). Indeed, using the method by Hladik et al. [17] we 
obtain 

A{A) = [3.7321, 6.7843] U [-0.6458, 0.3230] U [-4.1072, -1.0000]. 
Example 2. Consider the example given by Qiu et al. [27] (see also |18[ I35j l: 

/ [2975,3025] [-2015,-1985] \^ 

s_ [-2015,-1985] [4965,5035] [-3020,-2980] 

[-3020,-2980] [6955,7045] [-4025,-3975] 

\ [-4025,-3975] [8945,9055] / 

The local improvement (Algorithm [1]) yields an inner approximation 

/Xi(A^) = [12560.8377, 12720.2273], ti2{A^) = [7002.2828, 7126.8283], 
H^{A^) = [3337.0785, 3443.3127], fJ^^iA^) = [842.9251, 967.1082]. 

The vertex enumeration (Algorithmic]) produces the same result. Hence we can 
state that JIi{A^) and fi^{A^) are optimal. 

To call the last method, submatrix vertex enumeration (Algorithm |3]) we 
need an outer approximation. We use the following by [18] 

a;i(A^) = [12560.6296, 12720.2273], u}2{A^) = [6990.7616, 7138.1800], 
u}3{A^) = [3320.2863, 3459.4322], u;i{A^) = [837.0637, 973.1993]. 
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Now, submatrix vertex enumeration yields the same inner approximation as the 
previous methods. However, now we have more information. Since the outer 
approximation interval are mutually disjoint, the obtained results are the best 
possible. Therefore, ii^{A^) — Xi[A^), where i = 1, . . . ,4. 

Example 3. Herein, we present two examples for approximating the singular 
values of an interval matrix. Let A G ]^™x" ^nd q :— min{TO,n}. By the 
Jordan- Wielandt theorem ^ 151 EDI the singular values, <Ji{A) > • • • > aq{A)^ 
of A are identical with the q largest eigenvalues of the symmetric matrix 



A 

Thus, if we consider the singular value sets <Ti(A), . . . , (Tq{A) of some interval 
matrix A G M™^", we can identify them as the q largest eigenvalue sets of the 
symmetric interval matrix 

A"^' 



^■-\A 
(1) Consider the following interval matrix from [S]. 



A = 




Both the local improvement and the vertex enumeration result the same inner 
approximation, i.e. 

^ll{M) = [2.5616, 4.5431], /X2(M) = [1.2120,2.8541]. 

Thus, ai [A) ~ 4.5431. Additionally, consider the following outer approximation 
from dH]. 

uji{M) = [2.0489, 4.5431], U2{M) = [0.4239, 3.1817]. 

Using Algorithm [31 we obtain 

^l'l{M) = [2.5616, 4.5431], ^l'2{M) = [1.0000,2.8541]. 

Now we can claim that g_2{-^) = Ij since ujj^i^) > 0- Unfortunately, we cannot 
conclude about the exact values of the remaining quantities, since the two outer 
approximation intervals overlap. We know only that g_i{A) e [2.0489, 2.5616] 
and CT2(A) e [2.8541, 3.1817]. 

(2) The second example comes from Ahn & Chen ^2|. Let A be the following 
interval matrix 

/ [0.75, 2.25] [-0.015, -0.005] [1.7, 5.1] 
A = [3.55, 10.65] [-5.1, -1.7] [-1.95, -0.65] 
\ [1.05, 3.15] [0.005, 0.015] [-10.5, -3.5] 

Both local improvement and vertex enumeration yield the same result, i.e. 

^ll{M) = [4.6611, 13.9371], /X2(M) = [2.2140, 11.5077], 
/.t3(M) = [0.1296, 2.9117]. 
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Hence, cri(A) = 13.9371. As an outer approximation we use the following 
intervals, using |18) . 

a;i(M) = [4.3308, 14.0115], uj2{M) = [1.9305, 11.6111], 
uj3{M) = [0.0000, 5.1000]. 

Running the submatrix vertex enumeration, we get the inner approximation 

fi'iiM) = [4.5548, 13.9371], Mzl-^) = [2.2140, 11.5077], 
IJ-'siM) ^ [0.1296, 2.9517]. 

We cannot conclude that a^iA) — fJ-^{A) — 0.1296, because uj3{M) has a 
nonempty intersection with the fourth largest eigenvalue set, which is equal 
to zero. Also the other singular value sets remain uncertain, but within the 
computed inner and outer approximation. 

Notice that ^j/^{M) < ^^(M), whence fj/^{M) < Ai(M) = ^^(A) disproving 
the Hertz's Theorem 1 from [TB] that the lower and upper limits of Ai(Af) and 
A„(7Vf ) are computable by the vertex enumeration method. It is true only for 
Ai(M) and A„(M). 

Example 4. In this example we present some randomly generated examples of 
large dimensions. The entries of the midpoint matrix, Ac, are taken randomly in 
[—20, 20] using the uniform distribution. The entries of the radius matrix are 
taken randomly, using the uniform distribution in [0, i?], where i? is a positive 
real number. We applied our algorithm on the interval matrix M :— A^ A, 
because it has a convenient distribution of eigenvalue set — some are overlapping 
and some are not. Sharpness of results is measured using the quantity 

where e = (1,...,1)"^. This quantity lies always within the interval [0, 1]. The 
closer to zero it is, the tighter the approximation. In addition, if it is zero, then 
we achieved exact bounds. The initial outer approximation, u3i(M.^\ where 
1 < i < n, was computed using the method due of Hladik et al. [18J, and 
filtered by the method proposed by Hladik et al. in [T^- Finally, it was refined 
according to the comment in Section 14.3.21 For the submatrix vertex enumer- 
ation algorithm we implemented the branch & bound improvement, which is 
described in Section 14.3.11 and 14.3.21 
The results are displayed in Table [1] 

Example 5. In this example we present some numerical results on approxi- 
mating singular value sets as introduced in Example |31 The input consists of an 
interval (rectangular) matrix A C which is selected randomly as in the 

previous example. 

Table [2] presents our experiments. The time in the table corresponds to the 
computation of the approximation of only the q largest eigenvalue sets of the 
Jordan- Wielandt matrix. 
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n 


R 


Algorithm [T] 


Algorithm [2] 


Algorithm [3] 


sharpness 


time 


sharpness 


time 


sharpness 


time 


5 


0.001 


05817 


0.00 s 


05041 


0.00 s 


00000 


0.04 s 


5 


0.01 


0.07020 


0.00 s 


0.05163 


0.00 s 


0.00000 


0.03 s 


5 


0.1 


0.26273 


0.00 s 


0.23389 


0.00 s 


0.17332 


0.04 s 


5 


1 


0.25112 


0.00 s 


0.23644 


0.00 s 


0.20884 


0.01 s 


10 


0.001 


0.08077 


0.00 s 


0.07412 


0.09 s 


0.00000 


1.15 s 


10 


0.01 


0.13011 


0.01 s 


0.11982 


0.08 s 


0.04269 


1.29 s 


10 


0.1 


0.27378 


0.01 s 


0.25213 


0.09 s 


0.12756 


3.17 s 


10 


1 


0.56360 


0.01 s 


0.52330 


0.09 s 


0.52256 


2.58 s 


15 


0.001 


0.07991 


0.02 s 


0.07557 


7.3 s 


0.00000 


16.47 s 


15 


0.01 


0.21317 


0.02 s 


0.19625 


6.5 s 


0.11341 


2 m 29 s 


15 


0.1 


0.36410 


0.02 s 


0.34898 


7.0 s 


0.34869 


4 m 58 s 


15 


1 


0.76036 


0.02 s 


0.73182 


7.2 s 


0.73182 


7.5 s 


20 


0.001 


0.09399 


0.06 s 


0.09080 


7 m 21 s 


0.00000 


13 m 46 s 


20 


0.01 


0.24293 


0.06 s 


0.22976 


7 m 6 s 


0.12574 


1 h 14 m 55 s 


20 


0.1 


0.24293 


0.06 s 


0.22976 


7 m 14 s 


0.12574 


1 h 15 m 41 s 


20 


1 


0.82044 


0.06 s 


0.79967 


7 m 33 s 


0.79967 


7 m 39 s 


25 


0.001 


0.14173 


0.13 s 


0.13397 


6 h 53 m s 


0.02871 


9 h 32 m 54 s 



Table 1: Eigenvalues of random interval symmetric matrices A of dimension 
n X n. 



m 


n 


R 


Algorithm [T] 


Algorithm m 


Algorithm [3] 


sharpness 


time 


sharpness 


time 


sharpness 


time 


5 


5 


0.01 


0.08945 


0.00 s 


0.07716 


0.10 s 


0.00000 


0.53 s 


5 


5 


0.1 


0.09876 


0.01 s 


0.09270 


0.08 s 


0.00000 


0.73 s 


5 


5 


1 


0.43560 


0.01 s 


0.31419 


0.10 s 


0.26795 


4.34 s 


5 


10 


0.01 


0.11320 


0.02 s 


0.10337 


5.79 s 


0.00000 


7.91 s 


5 


10 


0.1 


0.13032 


0.02 s 


0.12321 


5.98 s 


0.00000 


8.40 s 


5 


10 


1 


0.35359 


0.02 s 


0.33176 


5.52 s 


0.22848 


21.53 s 


5 


15 


0.01 


0.10603 


0.05 s 


0.09424 


5 m 31 s 


0.00000 


5 m 36 s 


5 


15 


0.1 


0.17303 


0.04 s 


0.16758 


5 m 33 s 


0.00000 


7 m 58 s 


5 


15 


1 


0.46064 


0.05 s 


0.39708 


5 m 32 s 


0.31847 


15 m 47 s 


10 


10 


0.01 


0.10211 


0.06 s 


0.09652 


8 m 3 s 


0.00000 


8 m 19 s 


10 


10 


0.1 


0.13712 


0.07 s 


0.13387 


8 m 10 s 


0.00000 


14 m 12 s 


10 


10 


1 


0.39807 


0.07 s 


0.35580 


7 m 52 s 


0.30279 


26 h 48 m 38 s 


10 


15 


0.01 


0.09561 


0.12 s 


0.09116 


5 h 51 m 53 s 


0.00000 


5 h 54 m 56 s 



Table 2: Singular values of random interval matrices of dimension m x n. 
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6 Conclusion and future directions 

We proposed a new solution theorem for the symmetric interval eigenvalue prob- 
lem, which describes some of the boundary points of the eigenvalue set. Unfor- 
tunately the complete characterisation is still a challenging open problem. 

We developed an inner approximation algorithm (submatrix vertex enumer- 
ation), which in the case where the eigenvalue sets are disjoint, and the interme- 
diate gaps are wide enough, output exact results. To our knowledge, even under 
this assumption, this is the first algorithm that can guarantee exact bounds. 

Based on our numerical experiments suggest that the local search algorithm 
is superior to the submatrix vertex enumeration algorithm when the input ma- 
trices are not of very small dimension. 

Acknowledgment. ET is partially supported by an individual postdoctoral 
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